In [1]:
from IPython.core.display import HTML
css_file = '../ipython_notebook_styles/ngcmstyle.css'
HTML(open(css_file, "r").read())
Out[1]:
The aim of this computer laboratory is to gain a practical insight into the fundamental steps required to program a simple finite element code to solve a problem of elasticity. In doing so, you will learn:
The numerical environment chosen for this laboratory is IPython.
Let's consider a 2D solid body $\Omega $ at rest, subjected to body forces ${\bf{b}}$, prescribed displacements ${\bf{\bar d}}$ and tractions ${\bf{\bar t}}$. The problem of elasticity is defined as follows:
$$ \begin{align} {\text{div}}\sigma + {\bf{b}} & = {\bf{0}} & {\bf x}(x,y) &\text{ in } \Omega \\ {\bf{d}}(u,v) &= {\bf{\bar d}}(u,v) & {\bf x}(x,y) &\text{ on } \partial {\Omega _u} \\ {\bf{t}}({\bf x}) &= {\bf{\bar t}}({\bf x}) & {\bf x}(x,y) &\text{ on }\partial {\Omega_\sigma} \end{align} $$$\sigma $ is the symmetric Cauchy stress tensor.
Apply integration by parts to the PDE ${\text{div}}\sigma + {\bf{b}} = {\bf{0}}$ with the divergence theorem (relation between the flow (flux) of a vector field though a surface and the behavior inside the surface) recast it into its variational form. You will also take into account the conditions about the test function $\upsilon$ on the boundary where ${\bf{d}}$ is prescribed.
To relate stress and strain one must define what is called a constitutive relationship. In other branches of physics other types of relationships would be used (e.g. relation between strain and temperature).
In the case of isotropic linear elasticity (Hooke's elasticity) the Cauchy stress tensor is related to the strain tensor $\varepsilon $ through the following relationship:
$$ \begin{equation} \sigma = \mathbb{C} : {\bf \varepsilon} \end{equation} $$where $\mathbb{C}$ is the elasticity tensor.
Equation (2) is a tensor equation (valid in any coordinate system). If one decides to use a Cartesian coordinate system (2) can be written in index notation as:
$$ \begin{equation} {\left( {\bf \sigma} \right)_{ij}} = {\sigma _{ij}} = \left( \mathbb{C} : {\bf \varepsilon} \right)_{ij} = \mathbb{C}_{ijkl} \varepsilon_{kl} \quad \left\{ {{\text{2D}}:i,j,k,l:1 \dots 2; \quad {\text{3D}}:i,j,k,l:1 \dots 3} \right\} \end{equation} $$Please note that we use the following convention: bold characters represent vectors, matrices or tensors.
Because the Cauchy, strain and elasticity tensors are symmetric and the material isotropic one can simplify the tensorial constitutive equation into a matrix equation involving vectors and asecond-order matrix.
$$ \begin{equation} \left\{ {\bf \sigma} \right\} = \left\{ \begin{array}{l} {\sigma _{xx}}\\ {\sigma _{yy}}\\ {\sigma _{xy}} \end{array} \right\} = \underbrace {\frac{E}{{1 - {\nu ^2}}}\left[ {\begin{array}{*{20}{c}} 1&\nu &0\\ \nu &1&0\\ 0&0&{\frac{{1 - \nu }}{2}} \end{array}} \right]}_{\left[ \mathbb{C} \right]}\underbrace {\left\{ \begin{array}{l} {\varepsilon _{xx}}\\ {\varepsilon _{yy}}\\ 2{\varepsilon _{xy}} \end{array} \right\}}_{\left\{ {\bf \varepsilon} \right\}} \end{equation} $$$\left\{ {E,\nu } \right\}$ are material properties, respectively the Young's modulus and Poisson's ratio which characterise the stiffness and compressibility of the material.
Let’s consider a 4-noded bi-linear interpolation element (Figure 1).
The shape functions chosen here are simple linear Lagrange polynomials but other forms could be used (quadratic Lagrange, Hermite, NURBS, etc):
\begin{equation} {N_i} = \frac{1}{4}(1 + {\xi_i} \xi )(1 + {\eta _i}\eta ) \quad {\xi _i},{\eta _i} = - 1,1 \end{equation}There are as many shape functions as nodes in the element, i.e. 4:
\begin{equation} \left\{ \begin{array}{rl} N_1 &= \frac{1}{4}\left( {1 - \xi } \right)\left( {1 - \eta } \right)\\ N_2 &= \frac{1}{4}\left( {1 + \xi } \right)\left( {1 - \eta } \right)\\ N_3 &= \frac{1}{4}\left( {1 + \xi } \right)\left( {1 + \eta } \right)\\ N_4 &= \frac{1}{4}\left( {1 - \xi } \right)\left( {1 + \eta } \right) \end{array} \right\} \end{equation}We are using an isoparametric formulation which means that the position vector and displacement vector are interpolated using the same order of shape (or interpolation) functions.
where
$$ \begin{equation} {{\bf{x}}_A}(t) = {\rm{position}}\,(node\,A) = \left\{ \begin{array}{l} {x_A}\\ {y_A} \end{array} \right\} \end{equation} $$$$ \begin{equation} {\bf{x}}\left( {{\xi ^e},t} \right) = {\left\{ {\begin{array}{*{20}{c}} {x(t)}\\ {y(t)} \end{array}} \right\}_{node\,1}}{N_1}\left( {{\xi ^e}} \right) + {\left\{ {\begin{array}{*{20}{c}} {x(t)}\\ {y(t)} \end{array}} \right\}_{node\,2}}{N_2}\left( {{\xi ^e}} \right) + {\left\{ {\begin{array}{*{20}{c}} {x(t)}\\ {y(t)} \end{array}} \right\}_{node\,3}}{N_3}\left( {{\xi ^e}} \right) + {\left\{ {\begin{array}{*{20}{c}} {x(t)}\\ {y(t)} \end{array}} \right\}_{node\,4}}{N_4}\left( {{\xi ^e}} \right) \end{equation} $$where
$$ \begin{equation} {{\bf{d}}_A}(t) = {\rm{displacement}}\,(node\,A) = \left\{ \begin{array}{l} {u_A}\\ {v_A} \end{array} \right\} \end{equation} $$We can arrange all the degrees of freedom of the element in one vector:
$$ \begin{equation} {\bf{d}} = {\{ {u_1},{v_1},{u_2},{v_2},{u_3},{v_3},{u_4},{v_4}\} ^T} \end{equation} $$The discretised versions of the position and displacement vector can be rewritten in matrix form:
$$ \begin{equation} {\bf{x}}\left( {{\xi ^e},t} \right) = {N_I}\left( {{\xi ^e}} \right){{\bf{x}}_I}(t) = {{\bf{N}}^T}\left( {{\xi ^e}} \right){{\bf{x}}_I}(t) \end{equation} $$$$ \begin{equation} {\bf{d}}\left( {{\xi ^e},t} \right) = {N_I}\left( {{\xi ^e}} \right){{\bf{d}}_I}(t) = {{\bf{N}}^T}\left( {{\xi ^e}} \right){{\bf{d}}_I}(t) \end{equation} $$where
$$ \begin{equation} {{\bf{N}}^T}\left( {{\xi ^e}} \right) = \left\{ {{N_1}({\xi ^e}),{N_2}({\xi ^e}),{N_3}({\xi ^e}),{N_4}({\xi ^e})} \right\} \end{equation} $$The test function used for the formulation of the weak form is also discretised:
$$ \begin{equation} \upsilon \left( {{\xi ^e},t} \right) = {N_I}\left( {{\xi ^e}} \right){\upsilon _I}(t) = {{\bf{N}}^T}\left( {{\xi ^e}} \right){\upsilon _I}(t) \end{equation} $$Inject the discretised version of the displacement into the weak form established in Question 1. You will have to express the strain matrix using the ${\bf B}$ matrix (derivatives of shape functions). You should obtain 3 integral matrix/vector forms: the element stiffness matrix, the nodal load vector due to body forces and the nodal load vector due to boundary forces.
The integral expressions obtained in Question 2 will have to be numerically integrated using, in our case, Gauss integration. We have integrals defined over the element domain ${\Omega ^e}$ and its boundary $\partial {\Omega ^e}$.
There are well known integration formulas over parametric domains that are exact for polynomials like the Lagrange polynomials used in the shape functions. The integral of a function is transformed into the sum of weighted values that the function takes at specific points, the so-called integration points or Gauss points.
$$ \begin{equation} \int\limits_ {f(\xi )\,\,} d \xi = {w_i}f({\xi _i}) \end{equation} $$These integration formulas apply to the parent element (Figure 1) which is the element in the parametric domain $ = \{ \{ - 1,1\} ,\{ - 1,1\} \} $ so it is necessary to map back the parent element domain to the current (deformed) domain. This is a simple change of variables:
$$ \begin{equation} d{\Omega _e} = dx\,dy = \det ({\bf{J}})\,d\xi\, d\eta \end{equation} $$$$ \begin{equation} \int\limits_{{\Omega^{e}}} {f({\bf{x}})\,\,} d\Omega = \int\limits_ {f({\bf{x}}({\xi ^e}))\,\underbrace {\det \left( {\frac{{\partial {\bf{x}}}}{{\partial {\xi ^e}}}} \right)}_{{J_{{\xi ^e}}} = \det ({\bf{J}})}\,} d \xi = {w_i}{J_{{\xi ^e}}}({\xi ^e}_i)f({\xi ^e}_i) \end{equation} $$Here, we will use a 2 points Gauss integration rule: i.e. 4 integration points for the quadrilateral element:
$$ \begin{align} {{\bf{G}}_1} &= \left\{ { - \frac{1}{{\sqrt 3 }}, - \frac{1}{{\sqrt 3 }}} \right\} & {w_1} &= 1\\ {{\bf{G}}_2} &= \left\{ {\frac{1}{{\sqrt 3 }}, - \frac{1}{{\sqrt 3 }}} \right\} & {w_2} &= 1\\ {{\bf{G}}_3} &= \left\{ {\frac{1}{{\sqrt 3 }},\frac{1}{{\sqrt 3 }}} \right\} & {w_3} &= 1\\ {{\bf{G}}_4} &= \left\{ { - \frac{1}{{\sqrt 3 }},\frac{1}{{\sqrt 3 }}} \right\} & {w_4} &= 1 \end{align} $$Write a Python programme to:
Using the Python programme from question 3, solve the following problem.
Consider a 2D rectangular solid domain (length $L_x = 1m$; height $L_y = 0.5m$) with material properties $\{ E = 210{\rm{ GPa}};\,\,\,\nu = 0.3\} $ and sides labelled (counter-clockwise, starting from the left vertical side) as A, B, C and D.
The boundary conditions are as follows:
\begin{equation} u({\text{side}}\,{\text{A}}) = v({\text{side}}\,{\text{B}}) = 0 \end{equation}The load condition is:
\begin{equation} {\bf{F}}({\text{side C}})[Newton] = \left\{ \begin{array}{l} {F_x}({\text{side C}})\\ {F_y}({\text{side C}}) \end{array} \right\} = \left\{ \begin{array}{l} {10^6}\\ 0 \end{array} \right\} \end{equation}The boundary and loading conditions are homogeneous and therefore do not introduce shear coupling. The Cauchy stress tensor has only one non-null component, ${\sigma _{xx}}$.
NB: the total load of 10$^{6}$ N must be uniformly distributed over all the nodes of side C.
The strains can be calculated as
$$ \begin{align} {\varepsilon _{xx}} = {\sigma _{xx}}/E = {F_x}/E{L_y} &= - {\varepsilon _{yy}}/\nu \\ {\varepsilon _{xy}} &= 0 \end{align} $$and the displacements as
$$ \begin{align} u &= {\varepsilon _{xx}}{L_x} \\ v &= {\varepsilon _{yy}}{L_y}. \end{align} $$Do the results surprise you? Why?